; create a Surface density map from particle position (x,y,z) 
; to be used with sig2lens gravlens command

m=8.6e8
rhoc=2.7755e11

;if (n_elements(read) eq 1) then begin
;print, 'Data already loaded'
;endif else begin

;readcol,'halo388599.out',x1,y1
;readcol, 'halo1402911.out',x2,y2
;readcol, 'halo15.out',x3,y3

readcol, 'halo10.out',x4,y4,z4

; r200 (Mpc/h)
; 0.863868     0.764250     0.914033     0.951089

;read=1
;endelse

ind=where(sqrt(x4^2+y4^2) lt 0.2 and sqrt(x4^2+y4^2+z4^2) lt 0.95)

; Gravlens needs positions in kpc/h (see gravlens sig2lens command)

x=x4[ind]*1e3
y=y4[ind]*1e3

print, "N=",n_elements(x)

; set here the halo region to crop (with respect to the center position)

xmin=min(x)
xmax=max(x)
ymin=min(y)
ymax=max(y)

x=x+abs(xmin)
y=y+abs(ymin)

y=y/(ymax-ymin)
x=x/(xmax-xmin)

nbins=50
;binsize=(xmax-xmin)/float(nbins)
;print, "Crop region: dx=", xmax-xmin, " dy=", ymax-ymin, " binsize=", xbinsize, "x", ybinsize

;im = HIST_2D(x,y, MIN1=xmin, MAX1=xmax, MIN2=ymin, MAX2=ymax, BIN1=binsize, BIN2=binsize)

value=findgen(n_elements(x))
value[*]=m

im=tsc(value,x*nbins,nbins+1,y*nbins,nbins+1)

;use adaptative smoothing?

;result=asmooth_fft(im,2.5,plot=0)

; plot results

!p.charthick=2
!p.color=0
!p.background=-1

loadct,3

set_plot, 'ps'

;window,0,retain=2

device,/encapsulated,filename='sigma.eps',/inches,xsize=5.3,scale_factor=1.0,font_size=12,/bold,/color,bits=8

;tvscl,alog10(result.asm_im+0.05)

shade_surf,im
device,/close

; write data to file, this format can be read from gravlens sig2lens command

openw, 1, 'sigma.txt'
printf, 1,  nbins,' ', nbins
for i=1,nbins-1 do for j=1,nbins-1 do printf, 1, im[i,j]  
close,1

openw, 2, 'sigma.range'
printf, 2, xmin ,' ', xmax
printf, 2, ymin ,' ', ymax
close,2


end

